!(pwd() in LOAD_PATH) ? push!(LOAD_PATH, pwd()) : nothing
import LinearAlgebra
import Utilities

function expandedCheb(ngp)
    res = zeros(ngp)
    if ngp == 1
        res[1] = 0.0
    else
        for k = 1:ngp
            res[k] = cos((2 * (ngp - k + 1) - 1.0) * pi / (2 * ngp)) / cos(pi / (2 * ngp))
        end
    end
    return res
end

function dqcoeff!(xlist::Vector{Float64}, maxorder::Int64, c::Array{Float64,3})
    #xlist: input, one-dimensional real list, with length n
    #maxorder : input, interger
    #c : output, differential quadrature coefficient array, with dimensions (maxorder+1, n, n)
    n = size(xlist, 1)
    c .= 0.0
    for i = 1:n
        for j = 1:n
            if i != j
                c[1, i, j] = Mp(xlist, i) / ((xlist[i] - xlist[j]) * Mp(xlist, j))
            end
        end
    end

    for i = 1:n
        c[1, i, i] = -sum(c[1, i, 1:i-1]) - sum(c[1, i, i+1:n])
    end

    for m = 2:maxorder
        for i = 1:n
            for j = 1:n
                if i != j
                    c[m, i, j] = m * (c[m-1, i, i] * c[1, i, j] - c[m-1, i, j] / (xlist[i] - xlist[j]))
                end
            end
        end
        for i = 1:n
            c[m, i, i] = -sum(c[m, i, 1:i-1]) - sum(c[m, i, i+1:n])
        end
    end

    for i = 1:n
        c[maxorder+1, i, i] = 1.0
    end
    return nothing
end

function Mp(xlist, i)
    #implicit none
    #real(kind=iwp) :: Mp, xlist(:)
    #integer(4) :: n, i, j
    n = size(xlist, 1)
    r = 1.0
    for j = 1:(i-1)
        r = r * (xlist[i] - xlist[j])
    end
    for j = (i+1):n
        r = r * (xlist[i] - xlist[j])
    end
    return r
end

function LocalDQCoefficientsOdd(arr, L; maxorder=2)
    
    n = length(arr)
    ret = zeros(maxorder + 1, n, n)
    if L > n 
        println("Error! too large local point number L. L should be no more than total node number.")
        return ret
    elseif mod(L,2) == 0
        println("Error! L should be odd number.")
        return ret
    end

    c = zeros(maxorder+1, L, L)
    dqcoeff!(arr[1:L], maxorder, c)
    for m = 1 : maxorder + 1
        ret[m, 1 : (L + 1) ÷ 2, 1 : L] = c[m, 1 : (L + 1) ÷ 2, :]
    end
    for i = 1 + (L + 1) ÷ 2 : n - (L + 1) ÷ 2
        dqcoeff!(arr[i - (L - 1)÷2 : i + (L - 1)÷2], maxorder, c)
        for m = 1 : maxorder + 1
            ret[m, i, i - (L - 1)÷2 : i + (L - 1)÷2] = c[m, (L + 1)÷2, :]
        end
    end

    dqcoeff!(arr[n - L + 1 : n], maxorder, c)
    for m = 1 : maxorder + 1
        ret[m, 1 + n - (L + 1)÷2 : n, n - L + 1 : n, :] = c[m, (L + 1)÷2 : end, :]
    end
    return ret
end


function LocalDQCoefficientsEven(arr, L; maxorder=2)
    n = length(arr)
    ret = zeros(maxorder + 1, n, n)
    if L > n
        println("Error! too large local point number L. L should be no more than total node number.")
        return ret
    elseif mod(L, 2) == 1
        println("Error! L should be even number.")
        return ret
    end

    c = zeros(maxorder+1, L, L)
    dqcoeff!(arr[1:L], maxorder, c)
    for m = 1:maxorder+1
        ret[m, 1:1+L÷2, 1:L] = c[m, 1:1+L÷2, :]
    end
    for i = 2+L÷2:n-L÷2+1
        dqcoeff!(arr[i - L ÷ 2 : i + L ÷ 2 - 1], maxorder, c)
        for m = 1:maxorder+1
            ret[m, i, i-L÷2 : i + L ÷ 2 - 1] = c[m, 1+L÷2,:]
        end
    end

    dqcoeff!(arr[n-L+1:n], maxorder, c)
    for m = 1:maxorder+1
        ret[m, n-L÷2+2:n, n-L+1:n] = c[m, L÷2+2:end, :]
    end
    return ret
end

function LocalDQCoefficients(arr, L; maxorder=2)
    if mod(L,2)==0
        return LocalDQCoefficientsEven(arr, L, maxorder=maxorder)
    else
        return LocalDQCoefficientsOdd(arr, L, maxorder=maxorder)
    end
end

function test()
    n = 8
    maxorder = 3
    x = expandedCheb(n)
    a = 0.0
    b = 2.0 * pi
    x = 0.5 * (b-a) * x .+ (b+a) * 0.5
    c = zeros(maxorder + 1, n, n)

    # test GDQM coefficients
    dqcoeff!(x, maxorder, c)
    Utilities.saveMatrix(c[1,:,:], "Global_Order=1.txt", sep="\t")
    Utilities.saveMatrix(c[2,:,:], "Global_Order=2.txt", sep="\t")
    Utilities.saveMatrix(c[3,:,:], "Global_Order=3.txt", sep="\t")
    Utilities.saveMatrix(c[4,:,:], "Global_Order=0.txt", sep="\t")

    # test GDQM coefficients
    for L in 3:8
        c = LocalDQCoefficients(x, L, maxorder=3)
        Utilities.saveMatrix(c[1,:,:], "Local_Order=1"*"_bandWidth="*string(L)*".txt", sep="\t")
        Utilities.saveMatrix(c[2,:,:], "Local_Order=2"*"_bandWidth="*string(L)*".txt", sep="\t")
        Utilities.saveMatrix(c[3,:,:], "Local_Order=3"*"_bandWidth="*string(L)*".txt", sep="\t")
        Utilities.saveMatrix(c[4,:,:], "Local_Order=0"*"_bandWidth="*string(L)*".txt", sep="\t")
    end
end

test()